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Abstract. The study of diffusion and low frequency vibrational motions of particles 
on metal surfaces is of paramount importance; it provides valuable information 
on the nature of the adsorbate-substrate and the substrate-substrate interactions. 
In particular, the experimental broadening observed in the diffusive peak with 
increasing coverage is usually interpreted in terms of a dipole-dipole like interaction 
among adsorbates via extensive molecular dynamics calculations within the Langevin 
framework. Here we present an alternative way to interpret this broadening by means of 
a purely stochastic description, namely the interacting single adsorbate approximation, 
where two noise sources are considered: (1) a Gaussian white noise accounting for the 
surface friction and temperature, and (2) a white shot noise replacing the interaction 
potential between adsorbates. Standard Langevin numerical simulations for flat and 
corrugated surfaces (with a separable potential) illustrate the dynamics of Na atoms on 
a Cu(lOO) surface which fit fairly well to the analytical expressions issued from simple 
models (free particle and anharmonic oscillator) when the Gaussian approximation is 
assumed. A similar broadening is also expected for the frustrated translational mode 
peaks. 
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1. Introduction 

Diffusion and low-frequency vibrational motion of atomic and molecular adsorbates on 
surfaces are two very important and elementary dynamical processes. They provide 
valuable information on the nature of the adsorbate-substrate and adsorbate-adsorbate 
interactions. Moreover, a good understanding of these processes also constitutes a 
preliminary step in the study of more complicated phenomena in surface science, 
such as heterogeneous catalysis, crystal growth, lubrication, chemical vapor deposition, 
associative desorption, etc. Among the different experimental techniques used to study 
these processes and extract interaction potentials, [H El El HI [5] , the quasielastic He atom 
scattering (QHAS) [3l[6l[7l[8] is considered as the surface science analogue of quasielastic 
neutron scattering, which has been widely and successfully apphed to analyze diffusion 
in the bulk p[IO]. 

Diffusion times and thus diffusion coefficients may vary by orders of magnitude, 
depending on the temperature of the surface and the nature of the adsorbate-surface 
interaction (for a recent review, see for example, reference [1]). Scanning tunnelling 
microscopy (STM) and field ion microscopy (FIM) are especially useful for slow diffusion, 
where the time between jumps of the adatom is of the order of seconds. These 
measurements provide a direct observation of the diffusion process. QHAS is being 
used to determine diffusion processes in which the time between jumps is of the order 
of microseconds. In contrast with STM and FIM, QHAS provides information on 
the diffusion dynamics only indirectly, but it is one of the techniques where reliable 
adsorbate-adsorbate and adsorbate-substrate interactions can be obtained since the 
underlying theory is much more intrincate. Most of theoretical treatments applied to 
interpret QHAS results assume that the interaction of He atoms with the adsorbates 
and the surface is negligible. As far as we know, a general study on how the probe 
particles can distort the diffusion process is not available yet. 

In the QHAS technique, time-of-fiight measurements are converted to energy 
transfer spectra from which a wide energy range can be spanned and several peaks are 
observed. The prominent peak around the zero energy transfer, namely the quasielastic 
peak (Q peak), gives information about the adsorbate diffusion process. Additional 
weaker peaks at low energy transfers around the Q peak are also observed. These 
peaks are attributed to the parallel (to the surface) frustrated translational motion of 
the adsorbates (T modes) and also to the excitations of surface phonons (at positive 
energy transfers we have creation processes and annihilation ones at negative energy 
transfers). The measurable quantity experimentally is the so-called dynamic structure 
factor or scattering law, which gives the line shapes of all those elementary processes, in 
particular those corresponding to the Q and T peaks. The dynamic structure factor 
provides information about the dynamics and structure of the adsorbates through 
particle distribution functions, the latter being related to the nature of the adsorbate- 
substrate and the adsorbate-adsorbate interactions. 

At low coverages the interactions among adsorbates can be ignored, this allowing 
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to work within the so-called single adsorbate approximation. In this case, diffusion 
(or self- diffusion) is characterized by only taking into account the adsorbate-substrate 
interaction which is usually described by a stochastic model (Brownian-like motion) 
with two contributions: (1) the deterministic, phenomenological adiabatic potential, V , 
which models the interaction at T = 0; and (2) a stochastic force, Rcit), (Gaussian 
white noise) accounting for the vibrational effects induced by the temperature on the 
surface lattice atoms (and therefore on the adsorbates). In this framework, the dynamics 
is then carried out by solving the standard Langevin equation [H Ull |12l [13l [H] . In 
this type of dynamical calculation the parameters defining the interaction potential are 
chosen very often to obtain good agreement with the experimental data. This way of 
proceeding has the disadvantage that such potential functions are not unique, since there 
can be a multiplicity of forms leading to the same results; only ab initio calculations 
could provide, in principle, unique potential functions. 

When dealing with higher coverages, adsorbate-adsorbate interactions can no 
longer be neglected. In this case, the adsorbate-surface interaction can still be 
described as before, but pairwise potential functions accounting for the adsorbate- 
adsorbate interactions are usually introduced into Langevin molecular dynamics (LMD) 
simulations [8]. An alternative approach is to consider a purely stochastic description 
for these interactions, as we have shown elsewhere [151 IIS]- This description, which we 
have denominated the interacting single adsorbate approximation, is a (fully) Langevin 
approach based on the theory of spectral-line coUisional broadening developed by Van 
Vleck and Weisskopf [IT] and the elementary kinetic theory of gases [18]. This new 
approach explains fairly well experimental results where broadening is observed with 
increasing coverage in both the Q and the T-mode peaks [8]. The motion of a single 
adsorbate is modelled by a series of random pulses within a Markovian regime (i.e., 
pulses of relatively short duration in comparison with the system relaxation); the pulses 
simulate the collisions with other adsorbates. In particular, we describe the adsorbate- 
adsorbate collisions by means of a white shot noise as a limiting case of a colored shot 
noise [19]. The concept of shot noise has been applied to study thermal ratchets |20j . 
mean first passage times [21], or jump probabilities |22j; an interesting, recent account 
on general colored noise in dynamical systems can be found in |23j. 

In this work we present a detailed analytical and numerical description of the 
effects arising from both Gaussian white noise and shot noise when applied to surface 
diffusion problems. In particular, we show how the kinds of process taking place on 
a surface can be properly described analytically by means of a very simple theoretical 
framework, where the interplay between diffusion and low vibrational motions results 
very apparent. Accordingly, we have organized this work as follows. The description of 
the elements involved in the theoretical analysis of line shapes as well as the mathematics 
related to the analytical results are given in section [21 In section [3l we provide a 
numerical illustration of the surface dynamics for different types of models to give a 
broad perspective of our approach. Finally, the main conclusions arising from this work 
as well as some future trends are given in section HI 
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2. Stochastic approaches to surface diffusion 

2.1. Dynamic structure factor and intermediate scattering function 

In QHAS experiments one is usually interested in the differential reflection coefficient, 
which can be expressed as 



in analogy to scattering of slow neutrons by crystals and hquids [9l[T0]. This magnitude 
gives the probability that the probe He atoms scattered from the diffusing collective 
reach a certain solid angle VL with an energy exchange huj = Ej — Ei and a parallel 
(to the surface) momentum transfer AK = — Kj. In equation ([T]), is the 
(diffusing) surface concentration of adparticles; is the atomic form factor, which 
depends on the interaction potential between the probe atoms in the beam and the 
adparticles on the surface; and ^(AK,^;) is the dynamic structure factor or scattering 
law, which gives the line shapes of the Q and the T-mode peaks — other peaks can 
also be present, such as the inelastic ones related to surface phonon excitations — and 
provides a complete information about the dynamics and structure of the adsorbates 
through particle distribution functions [HI [10]. Experimental information about long 
distance correlations is obtained from the scattering law when considering small values 
of AK, while information on long time correlations is available at small energy transfers 

Surface diffusion is a dynamical problem that can be tackled by means of classical 
mechanics for heavy adparticles. Thus, let us consider an ensemble of classical 
interacting particles on a surface. Their distribution functions are described by means 
of the so-called van Hove or time-dependent pair correlation function G(R,t) pj. This 
function is related to the dynamic structure factor as 



Given an adparticle at the origin at some arbitrary initial time, G(R,t) represents 
the average probability for finding a particle (the same or another one) at the surface 
position R = [x, y) at a time t. Note that this function is a generalization of the 
well known pair distribution function g(R) from statistical mechanics flEl [2l] . since it 
provides information about the interacting particle dynamics. 

Depending on whether correlations of an adparticle with itself or with another one 
are considered, a distinction can be made between the se// correlation function, (7^(1^; t), 
and the distinct correlation function, Gd(R.,t). The full pair correlation function can 
then be expressed as 



According to its definition, Gs(R.,t) is peaked at t = and approaches zero as 
time increases because the adparticle loses correlation with itself. On the other 
hand, Gd(R.,t) = g{lV) gives the static pair correlation function — the standard pair 




(1) 




(2) 




(3) 
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distribution function — at t = 0, and approaches the mean surface number density cr of 
diffusing particles as t ^ oo. Accordingly, equation ([3]) can be split up as 

G(R,0) = 5(R)+^(R) (4) 

at t = 0, and expressed as 

G{R,t)^a (5) 

for a homogeneous system with ||R|| — »• oo and/or t — >■ oo. At low adparticle 
concentrations {9 <^ 1), i.e., when interactions among adsorbates can be neglected 
because they are far apart from each other, the main contribution to equation ([3]) is 
Gs (particle-particle correlations are negligible and ~ 0). On the contrary, for high 
coverages it is expected that Gd presents a significant contribution to equation ([3]). 

Within this theoretical framework, the dynamic structure factor is better written 
as [9] 



5(AK,w) = / e-*^*/(AK,t)rft, (6) 



where 



1 ^ 

/(AK,t) = i_(^e-^^^-^^(°^e*^^-^^'(*)) (7) 

hi' 

is the intermediate scattering function — note that this function is the space Fourier 
transform of (^(R, t). In equation (I7j) the brackets denote an ensemble average and 
Rj(t) the trajectory of the jth adparticle on the surface. This function can again be 
split up into two sums, distinct {Id) and self (/g), depending on whether the crossing 
terms are taken into account or not, respectively. Following the language used in neutron 
scattering theory, the corresponding Fourier transforms of / and Is give what is called the 
coherent scattering law ^(AK, uj) and incoherent scattering law 5's(AK, a;), respectively. 
In QHAS experiments, and with interacting adsorbates, coherent scattering is always 
obtained. The corresponding theoretical interpretation of that scattering is usually 
carried out in terms of the Vineyard convolution approximation [10] , where the distinct 
correlation function Gd is expressed as a convolution of the self correlation function Gg- 
This approximation is known to fail at small distances where the surface lattice becomes 
important. Whereas in neutron scattering many attempts to improve the convolution 
approximation have been developed, within the QHAS context very little effort has been 
devoted to this goal. 

At finite coverages, one usually distinguishes between two diffusion coefficients: the 
tracer diffusion constant {Dt) and the collective diffusion constant {Dc) [1]. The first one 
refers to the self-diffusion process and focuses on the motion of a single adsorbate. On 
the contrary, Dc is related to the collective motion of all adsorbates which is governed 
by Fick's law. In any case, a Kubo-Green formula relates Dt or D^ with the velocity 
autocorrelation function of a single adsorbate or with the corresponding for the velocity 
of the center of mass, respectively. 
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2.1.1. The dipole-dipole interaction For interacting adsorbates, LMD simulations 
are widely carried out. Typically, one solves numerically a system of coupled 
Langevin equations where denotes the number of adsorbates. In most of systems, 
the Markovian-Langevin approximation is assumed because the Debye energy of the 
substrate excitations is greater than the adsorbate T-mode energy and therefore the 
damping can be considered as instantaneous (memory effects are negligible). The 
adsorbate-adsorbate interaction is also assumed pairwise and is given by a repulsive 
dipole-dipole potential. The so-called Topping's depolarization formula relates the 
dipole strength with the coverage [25]. This type of interaction is attributed to the 
electrostatic repulsion between the dipoles due to the charge transfer from the adatoms 
to the substrate. Numerical discrepancies between the experimental and molecular 
dynamics results are found for the broadening of the Q and the T-mode peaks as a 
function of the parallel wave vector momentum transfer. The discussion about the 
origin of such discrepancies is still an open problem [8]. 

2.2. The interacting single adsorbate approximation 

When the coverage increases, the introduction of pairwise potential functions results in 
a realistic description of the adsorbate dynamics. However, there is no a simple manner 
to handle the resulting calculations by means of a simple theoretical model, and one 
can only proceed by using some suggested fitting functions. Moreover, carrying out 
LMD simulations always result in a relatively high computational cost due to the time 
spent by the codes in the evaluation of the forces among particles. This problem gets 
worse when working with long-range interactions, since a priori they imply that one 
should consider a relatively large number of particles in order to get a good numerical 
simulation. One way to avoid such inconveniences (interpretative and computational) 
could consist in using a simple, realistic stochastic model. For a good simulation of a 
diffusion process, one has to consider very long times in comparison to the timescales 
associated with the friction caused by the surface or to the typical vibrational frequencies 
observed when the adsorbates keep moving inside a surface well. This means that 
there will be a considerably large number of collisions during the time elapsed in the 
propagation, and therefore that, at some point, the past history of the adsorbate could 
be irrelevant regarding the properties we are interested in. This memory loss is a 
signature of a Markovian dynamical regime, where adsorbates have reached what we 
call the statistical limit. Otherwise, for timescales relatively short, the interaction is not 
Markovian and it is very important to take into account the effects of the interactions 
on the particle and its dynamics (memory effects). The diffusion of a single adsorbate 
is thus modeled by a series of random pulses within a Markovian regime (i.e., pulses of 
relatively short duration in comparison with the system relaxation) simulating collisions 
between adsorbates. In particular, we describe these adsorbate-adsorbate collisions by 
means of a white shot noise [19] as a limiting case of colored shot noise. In this way, 
a typical LMD problem involving adsorbates is substituted by the dynamics of a 
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single adsorbate, where the action of the remaining — 1 adparticles is replaced by a 
random force described by the white shot noise. In this approximation, the distinction 
between self and distinct time-dependent pair correlation function does not exist and 
equations ([2]) and ([6]) still hold. The intermediate scattering function reads now as 

/(AK, t) = (e-*^K.[R{t)-R(0)]^ ^ ^g-iAK./X*') dt'^ (g) 

In order to get some analytical results and therefore a guide for the interpretation of the 
numerical Langevin simulations, the intermediate scattering function can be expressed 
as a second order cumulant expansion in AK, 

/(AK, t) ^ e-^^'/oC^-OCAKCi'K'^ (9) 



where 



CAK(r) = (i;AK(0)t;AK(r)) = lim 1= [ v^^{t) v^^{t + t) dt (10) 

is the autocorrelation function of the velocity projected onto the direction of the parallel 
momentum transfer (whose length is AK = ||AK||). Only differences between two 
times, r, are considered because this function is stationary. This is the so-called 
Gaussian approximation [IB], which is exact when the velocity correlations at more 
than two different times are negligible, thus allowing to replace the average acting over 
the exponential function by an average acting over its argument. This approximation 
results of much help in the interpretation of numerical simulations as well as in getting 
an insight into the underlying dynamics. 

The decay of Cak(t) allows to define a characteristic time, the correlation time, 



oo 



r = -4r / CAK(T)(ir, (11) 



where \/Jv^ = a/ ksT /m is the average thermal velocity in one dimension — though 
the dimensionality is two, note that Cak is defined along a particular direction (that 
given by AK) and therefore Cak(O) = /csT/m, i.e., the square of the thermal velocity in 
one direction. The correlation time is related to the line shape broadening in the sense 
that it provides a timescale for the decay of the intermediate scattering function and 
therefore information about the width of the Q peak in the dynamic structure factor. 



2.2.1. Some elementary notions on Gaussian white noise Alternatively to the Einstein- 
Wiener stochastic model for Brownian motion, in 1930 Orstein and Uhlenbeck 
formulated another one taking the particle velocity as the stochastic variable of interest 
(rather than its position, which is found integrating the velocity). The starting point 
for this model is the equation proposed by Langevin in 1908, 

mv = —m'-fv + mRcit). (12) 

This the simplest expression for an equation describing the Brownian motion of a particle 
of mass m in one dimension. The right hand side (r.h.s.) of this equation can be split 
up into two contributions: (1) a deterministic part characterized by the friction force 
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—m'-fv (7 being the friction coefficient depending on the fluid viscosity) and (2) a random 
part governed by the random force mRcit), where Rq is a Gaussian white noise. The 
random force (or, equivalently, the Gaussian white noise) satisfies two conditions: 

(i) The process Rcif) is Gaussian with zero mean, i.e., {Rait)) = 0. 

(ii) The force-force correlation time is infinitely short, i.e., m? {Rg{^)Rg{t)) = ^^{'t)i 
A being a constant that gives the strength of the couphng between the adparticle 
and the environment. 

The validity of this model relies on the fact that the Brownian particle is much heavier 
than the particles constituting the environment. This implies that the kicks received 
by the particle of interest, though relatively weak, are very effective when considered 
in a very large number (the central limit theorem holds and therefore the noise is 
Gaussian). When applied to the motion of single adsorbates, the kicks come from 
the surface fluctuations at a given temperature. Note that the time evolution of the 
environmental degrees of freedom is not taken into account because their correlations 
decay faster than those of the particle (Markovian approximation), as expressed by the 
condition (ii). 

The relationship between the friction in the Langevin equation and the fluctuations 
of the random force is given by the fluctuation- dissipation theorem [26], which reads as 

7(0;) = / {SRciO) SRcir)) e— dr (13) 
ksT Jo 

where 

SRcit) = Rcit) - {Rait)) (14) 

is the fluctuation due to the random noise function Rcit). Making use of the properties 
(i) and (ii) for the Gaussian white noise, the r.h.s. of equation (fT3l) becomes 

That is, the frequency spectrum of the friction force is flat, or white in the sense that all 
frequencies contribute equally to this spectrum, in analogy to white light. This allows 
to establish 

7(^)=7, (16) 

with the strength of the coupling between the Brownian particle and the environment 
thus being 

A = 2mjkBT, (17) 
The Gaussian white noise correlation function can now be defined as 

Gcir) ^ {6Rg{0)6Rg{t)) = ^^5(r). (18) 

m 

This dynamics implies that at thermal equilibrium the equipartition theorem holds. 
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2.2.2. Some elementary notions on shot noise The concept of noise arises from the 
early days of radio, the so-called shot noise [27] being one of the main sources of noise, 
which was first considered by Schottky [28] in 1918. Studies on this type of noise were 
developed during the 1920 's and 1930 's, and were summarized and largely completed 
by Rice [29] in the mid 1940's. The paradigm of shot noise is the non steady electrical 
current generated by independent (i.e., non correlated) electrons arriving at the anode 
of a vacuum tube. This (time-dependent) electric current can be expressed as 

I{t) = Y,Ht-U), (19) 

i 

where the pulse function hi{t — ti) represents the contribution to the current due to the 
ith individual electron, which is assumed to be identical for each electron. Moreover, 
it is also assumed that each electron arrives independently of the previous ones. The 
arrival times ti are thus randomly distributed with a certain average number per unit 
time according to a Poisson distribution. 

In our case, the current of electrons is replaced by the impacts received by an 
adsorbate from other surrounding adsorbates. We express this random force as mSRs{t), 
where 

6Rs{t) = Rs{t) - {{Rs)) (20) 
after making use of equation iHM . and where 

= (21) 

K 

As seen, the double average bracket in the last expression indicates averaging over the 
number of collisions (K) according to a certain distribution {Pr) and the total time 
considered (T). In analogy to equation (fT9|) we have 

K 

Rs{t) = J2bk{t-tk). (22) 

k=l 

Here, bk(t — tk) provides information about the shape and effective duration of the kth 
adsorbate-adsorbate collision at tj.. Then the probability to observe K collisions after a 
time T follows a Poisson distribution [19] , 

PKiT) = ^^e-^^, (23) 

where A is the average number of collisions per unit time. That is, now the friction 
coefficient has to be interpreted in terms of the collision frequency between adsorbates. 

Assuming sudden adsorbate-adsorbate collisions (i.e., strong but elastic collisions) 
and that after-collision effects relax exponentially at a constant rate A', the pulses in 
equation (l22i) can be modeled as 

&fc(t-4)=CfcA'e-"'(*-*^), (24) 

with t — tk > and Ck giving the intensity of the collision impact. Within a realistic 
model, collisions take place randomly at different orientations and energies. Hence it is 
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reasonable to assume that the coefficients are distributed according to an exponential 
law, 

^(cfc) = -e-'='=/", c, >0, (25) 
a 

where the value of a will be determined later on. Independently of their intensity, 
within this model any pulse decays at the same rate A', as seen in equation fl2^ . This 
rate defines a decay timescale for collision events, = 1/A'. On the other hand, 
the (collision) friction coefficient introduces a new timescale = 1/A, which can be 
interpreted as the (average) time between two successive collisions. The diffusion motion 
of the interacting adsorbates will take a time of the order of r,. in getting damped. 

As previously done with the Gaussian white noise, we can now compute the time 
correlation function of the shot noise 

Gs{r) = {{5Rs{Q)5Rs{T))), (26) 

where the double bracket is defined as in equation fl2T]) . A general expression for Qsi.'^) 
can be readily obtained after straightforward algebraic manipulations, which yield 

Gs{r) = l-y^PK{r)K r {h{t-t')h{t + T-t')),dt\ (27) 
^ K 

where 

{■),= / ■gic)dc (28) 



Jo 

is the average over the pulse intensity, with g{c) given by equation (125|) . Taking into 
account that 'YIik ^k{T)K = AT and introducing the change of variable ( = t — t', 
equation fl271) can be written in an approximated way as 

/oo 
{biOKC + r))cdC. (29) 
-oo 

This approximation relies on the hypothesis that b{() ^ outside a narrow time interval, 
< C < A, with A a few times larger than Tc, but of the same order of magnitude. Thus, 
equation fl29|) is a general expression independent of the shape considered to simulate 
the pulses. Substituting equation into fl2I?]) one then obtains 

^5(r) = ^e-^'M. (30) 

The validity of the standard Langevin approach to study a dynamics governed by 
a shot noise is determined by the fluctuation-dissipation theorem [26] . As will be seen, 
in the generalized Langevin formulation we have a memory function in terms of the 
time-correlation function instead of a (time-independent) friction coefficient. According 
to the fluctuation-dissipation theorem, the formal relationship between the frequency 
spectrum of the memory function and the random force correlation function is given by 

POO 

^~(^) = r^ / Ss{r)e-^--dr. (31) 
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Introducing equation (15UI) into this integral yields 
whose real part is 



Mi('^)] = \ U'^) + C(^)] = A ^s^. (33) 
Two limits are interesting in this expression: \' <^ uo and A' ^ uo. The first limit 
involves very short timescales (smaller than Tc). Memory effects are important and the 
generalized Langevin equation should be applied. Note that in this case, equation ([32D 
can be written as 

fM^A^, (34) 

and this frequency-dependent friction does not allow to define an appropriate relaxation 
timescale (colored shot noise). Conversely, in the second case, the collision timescale 
rules the system dynamics; it establishes a cutoff frequency, which leads to 

e» ^ A f 1 - (35) 



This expression can be approximate as ^{oj) ~ A whenever \ <^ uj <^ ujc = (i.e., 
Tr <C Tc). This limit holds for strong but localized (or instantaneous) collisions (as 
assumed here) as well as for weak but continuous kicks (Brownian motion). Since it 
is similar to the condition leading to equation (fT6|) . we can speak about a Poissonian 
white noise and make use of the standard Langevin equation. 

Due to their influence in the diffusion process, it is worth comparing the time series 
corresponding to SRcit) and 6Rs{t). For the sake of simplicity, we consider that all 
pulses in the shot noise have the same intensity — it can be easily shown that this is 
equivalent to assume the pulse intensities distributed according to equation (!25l) . though 
substituting (cfc) and (c^) by C and (C = 2kBT/m being the impact strength), 
respectively, wherever such averages appear. In figure [11(a) we observe the time series 
corresponding to a Gaussian white noise, SRcit). This series presents a fractal-like 
profile, i.e., enlargements of any subinterval will look like the whole interval considered. 
This is because within any relatively short period of time the particle feels many 
(uncorrelated) kicks from the surface in both directions, positive and negative. This 
is not the case, however, for a shot noise with A = 7. As seen in figure [Uj^b), the particle 
receives 2 kicks in average from another adsorbates in a much larger period of time 
(about 50 times the one considered for SRcit)). This means that in principle the mean 
free path for a single adparticle is relatively large. However, though a single adparticle 
receives very few hits, the collective effect is similar to that of having a Gaussian white 
noise because of the random distributions of the kicks. The effect is more apparent 
when larger values of A are considered, because the mean time between consecutive 
collisions decreases dramatically, as can be seen in figure [U^c) . In any case, both noises 
are completely uncorrelated, that is, 

{SRG{t)SRs{t')} = 0. (36) 
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Figure 1. (a) Random noise function, SRcit), corresponding to a Gaussian "white 
noise -with 7 = 0.1 ('^o = 2.2049 x lO"** a.u. is the harmonic frequency associated 
"with the nonseparable adsorbate-substrate interaction potential (sec section [331 )• (b) 
Random noise function, 5Rs{t), corresponding to a shot noise "with A = 7. (c) Random 
noise function, SRs{t), corresponding to a shot noise "with A = 10 7. The calculations 
refer to Na atoms at T = 200 K and assuming A' = 10^^^ a.u. (41.3 ps ^^). 



2.2.3. Relationship between coverage and the collision frequency Here we "would like to 
show ho"w the coverage 9 and A can be related in a simple manner. In the elementary 
kinetic theory of transport in gases [18] diffusion is proportional to the mean free path 
Z, "which is proportionally inverse to both the density of gas particles and the effective 
area of collision "when a hard-sphere model is assumed. For two-dimensional collisions 
the effective area is replaced by an effective length (twice the radius p of the adparticle) 
and the gas density by the surface density a. Accordingly, the mean free path is given 
by 

1 = —^. (37) 

Taking into account the Chapman- Enskog theory for hard spheres, the self-diffusion 
coefficient can be written as 



1 keT , , 

= — a/^. 38 
opo" V "in 
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Now, from Einstein's relation (see section [2.2.4p . and taking into account that 9 = a? a 
for a square surface lattice of unit cell length a, we finally obtain 



Therefore, given a certain surface coverage and temperature, A can be readily estimated 
from equation (l39l) . Notice that when the coverage is increased by one order of 
magnitude, the same holds for A at a given temperature. 

2.2.4- Interacting adsorhate dynamics within the Langevin equation Except for the 
case of a free-potential regime, the analytical study of particle motion in two dimensions 
results intractable in general due to the correlations between the two degrees of freedom. 
However, since we are interested in offering an analytical formulation that allows to 
better understand the process ruling surface diffusion and low vibrational motions, it is 
sufficient to proceed in one dimension and then try to adapt the resulting formulation 
to two dimensions. Thus, the motion of an adsorbate subjected to the action of a bath 
consisting of another adsorbates on a static one-dimensional surface potential can be 
well described by a generalized Langevin equation, 

x{t) = - [ 7]{t- t') x{t') dt' + F{x{t)) + 6R{t), (40) 
Jo 

where x represents the adsorbate coordinate, and its first and second time derivatives 
are expressed by one and two dots on the x variable, respectively. In this equation, 7]{t) 
is the bath memory function, which includes the effects arising from both the Gaussian 
white noise and the shot noise. Thus, the noise source is expressed as the sum of both 
contributions, 

6RGs{t) = 6RG{t) + 6Rs{t). (41) 

Regarding the deterministic term, F = —W is the deterministic force per mass unit 
derived from the periodic surface interaction potential iV{x) = V{x + a), a being the 
period along the x direction). If Tc is relatively small (i.e., collision effects relax relatively 
fast), the memory function in equation ( HOl) will be local in time. This allows to express 
the memory function asri{t — t') ~ (7 + A)5(t — t') and expand the upper time limit in the 
integral to infinity. Taking into account this approximation and from the fiuctuation- 
dissipation theorem, which leads to r] = •y + X, equation (HOl) becomes 

x{t) = -r]x{t) + F{x{t)) + dRcit) + 6Rs{t). (42) 

The solutions of this equation can be straightforwardly obtained by formal integration, 
this rendering 

v{t) = voc-"" + f e-''(*-*')F(x(t')) dt' + f e-"'^'-''^ 5 Resit') dt', (43) 
Jo Jo 

x{t) = xo + - (1 - e-"*) + - / [l - e-''^*-*')] F{xit')) dt' 
V V Jo ^ J 



1 '■^ 



+ - 
V Jo 



1 - 



6RGsit')dt', (44) 
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where Vq = f (0) and Xq = x(0). As can be seen, for SRgs = 0, equations and 
become the formal solutions of purely deterministic equations of motion. Therefore, 
without loss of generality, these solutions can be more conveniently expressed as 

V{t) = Vd{t) + Vs,G{t) + Vs,s{t), (45) 
X{t) = Xd{t) + Xs,G{t) + Xs,s{t), (46) 

where d refers to the deterministic terms and s to those depending on the stochastic 
forces. Nonetheless, note that when 6 Resit) 7^ the deterministic part will also present 
some stochastic features due to the evaluation of F{x) along the trajectory x(t), which 
is stochastic. 

Taking advantage of equation ( l36l) . we can write 

ivit)) =Vdit), (47) 

{v\t))=vm + {vlait)) + (visit)), (48) 
(x(t)) =xait), (49) 
(x'it)) = xlit) + {xlait)) + (xlsit)), (50) 

where the barred magnitudes indicate the respective averages of the deterministic part 
of the solutions, and 



Kxit)) 



-2r]t 



t-t' 



1 



Kxit)) = - , 

T Jo 



dt' 



-v{t-t') 



e'^Gxir) dr, 
t~t' 



^■q{t-t'-T) 



(51) 



Gxir)dr, (52) 



with X = G oi S. For X = G, the final form of these expressions reads as 



Kcit)) 

Kait)) 
and for X = S* as 

(visit)) 

«sit)) 



2l 

7f 



1-e 



-27)1 



2r]t+ 1 



-fit 



A A' 



Vrj (A' — 7]) 



+ 



A 1 

2 

2rf 



Irit 



4A' 



1-e" 
- 2ri 



-277* 



A'-r/ 



2A^A 

' A'2 - 



1-e 
2A' 



-(A'+r?)t 



(53) 
(54) 

(55) 



2rf 



A'(A' - T]) 

Introducing now the assumption A' ^ A, we obtain 

A 1 



A'2 - if 



2iX'-v) 



1-e 



-2rit 



-2r]t 



2r]t+l 



2 - e"' 



(56) 

(57) 
(58) 



These equations are identical to equations (!53|) and (!54|) . except for a weighting factor 
X/t]. The weighting factors 'y/rj and X/t] indicate the contribution arising from each 
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noise. Adding equations (1^ and (157|1 . on the one hand, and equation (IMI) and 
on the other hand, we reach 

{vl(t))^^^[l-e-'^'), (59) 

(^sW)-^^[2r^t+l-(2-e-''*)'], (60) 

where no clue about the type of noise is left. Equation (|59ll can now be used to determine 
the value of a; assuming that the equipartition theorem satisfies for t ^ cxd, 

\m{v\^)) = ]^kBT. (61) 

Therefore, since v'j^{t) = v^e"'^'^^ and the time-dependent term in equation (l59l) vanishes 
asymptotically. 



"^V^il- ''''' 

If the system is initially thermalized (i.e., it follows a Maxwell-Boltzmann distribution 
in velocities) and has a uniform probability distribution in positions around x = 0, then 
vq = 0, -Uq = kBT/m and xq = 0. Thus, for A' ^ A (i.e., in the limit of a Poissonian 
white noise), we finally obtain 

Ht)) =0, (63) 
(At)) = (64) 



m 



{x{t)) =0, (65) 



knT 



{x^it)) = xl + ^ 2'nt+l-(2-e-^') . (66) 
m?7^ L V / J 

Obviously, these equations constitute a limit and therefore for values of the parameters 
out of the range of the approximation deviations will be apparent. 

As happens with Brownian motion, two regimes are clearly distinguishable from 
equation ( 166|) . For rjt collision events are rare and the adparticle shows an almost 
free motion with relatively long mean free paths. This is the ballistic or free- diffusion 
regime, characterized by 

{x'{t))^^-^t\ (67) 
m 

On the other hand, for rjt ^ 1 there is no free diffusion since the effects of the 
stochastic force (collisions) are dominant. This is the diffusive regime, where mean 
square displacements are linear with time: 

^^2^t))r^^^t = 2Dt. (68) 
mrj 

This is the so-called Einstein's law. Note from equation fl68l) that: (1) lowering the 
friction rj acting on the adparticle leads to a faster diffusion (the diffusion coefficient D 
increases) and (2) diffusion becomes more active when the surface temperature increases. 
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Figure 2. Upper panels: Trajectory dynamics for V{x, y) = and: (a) a Gaussian 
wlrite noise with 7 = 0.1 (wq — 2.2049 x 10~^ a.u.) and (a') a shot noise with 
A = 7. Lower panels: Phase diagram corresponding to the trajectories represented 
above. The trajectories have been propagated for 2.5x10^ a.u. and A' = 10^"^ a.u. 



To conclude this section, it is worth comparing the type of dynamics ruled by each 
noise, Gaussian or shot. In figures [2] and [3] we present some trajectories together with 
their phase diagram (for the x-coordinate) for V{x, y) = and V{x, y) as given by 
equation ( 198|) (see section [3]), respectively. The noise functions are such that A = 7 
(figures [T](a) and [U^b)) and the numerical details are the same as in figure [1] (the 
total propagation time is t = 2.5 x 10^ a.u.). The first striking difference that we 
can appreciate is that, because of the continuous kicking in the case of a Gaussian 
white noise, the trajectory looks much smoother in the case of a shot noise. This can 
be better seen in the upper panels of figure [3] — the effect is not so pronounced in 
figure [2] because we have only represented one out of 50 time steps — and also in the 
lower panels of both figures. The lack of smoothness in the trajectories affected by the 
Gaussian white noise is a manifestation of this type of kicking (see figure [T](a)). On 
the other hand, it is also interesting to observe the difference in the dynamics induced 
in both cases when a potential is introduced. The presence of potential wells, where 
a particle can remain trapped, will make the particle to undergo a vibrational motion 
inside such wells. The low frequency vibrational motion observed gives rise to the so- 
called frustrated translational mode or T mode. 

2.2.5. The velocity autocorrelation function The average magnitudes given in 
section [2.2.41 are only relevant to understand the ensemble behaviour of the trajectories 
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Figure 3. Upper panels: Trajectory dynamics for V{x, y) as given by model 3 (see 
section [SH]) and: (a) a Gaussian white noise with 7 — 0.1 and (a') a shot noise with 
A = 7. Lower panels: Phase diagram corresponding to the trajectories represented 
above. The trajectories have been propagated for 2.5x10^ a.u. and A' = 10^"^ a.u. 



along time, and therefore to have some idea on the dynamics induced by the 
corresponding interaction potential. In order to get a deeper insight into the diffusion 
process — and therefore also into the line shapes associated with the dynamic structure 
factor — it is important first to compute the velocity autocorrelation function given by 
equation (ITOl) . This function can be expressed as 

ft ft+T — t' 

C{t) = {v^^{t)v^^{t + T)) = e-''^-^^ \ dt'e^^'' / e'^'GGs{s)ds.{m) 

Jo J~t' 

As seen, this expression only depends on the stochastic part of the solution for the 
velocity; the deterministic one cancels out in the averaging procedure. Moreover, it 
also takes into account that there are no correlations among the two space coordinates 
[x and y). This is true in any case for the Gaussian white noise; for the shot noise it 
only holds in the white noise limit, otherwise some correlations among the two degrees 
of freedom could still be present when V due to the finite duration of the pulses. 
To avoid this inconvenience, which deviates the treatment from analyticity, here we 
assume conditions for which the white noise limit is satisfied (see section [2. 2. 4p . Finally, 
note that in equation fl69|) it is also considered that the motion along both coordinates 
contribute equally (i.e., {vx{0) Vx{t)) ~ {'VyiO) Vy(t))). This finds its justification in the 
fact that the potential models that will be used display the same symmetry along each 
direction and therefore the corresponding dynamics is expected to be the same (this 
would not be the case, for example, when the activation barrier is much higher in one 
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direction). In this sense, the calculation of the velocity autocorrelation function can be 
regarded as one- dimensional. 

The analytical results derived following the assumptions described above will 
be compared with the numerical simulations, which have been obtained using the 
corresponding exact expressions (e.g., the r.h.s. of the first equality in equation ([H])). In 
order to better interpret our exact numerical results two limiting cases will be studied 
below. It is precisely the discussion in terms of both cases what will allow us to better 
understand the physical processes taking place on the surface. 



2.2.6. The low corrugation regime. Quasi-free adparticles. In the case of diffusion on 
low corrugated surfaces, where the role of the adiabatic adsorbate-substrate interaction 
potential is negligible and only the action of the thermal phonons is relevant, one can 
assume V ^ 0. The adparticle motion can then be regarded as quasi- free since it is not 
ruled by a potential, but only influenced by the stochastic forces. Within this regime, 
equation (!69l) becomes 

Cir) = Itsle-"^ + ^^^T^ (e- - ^ e--) . (70) 
rj m rj m A'^ — rj-^ \ a / 

In this expression the contributions arising from each noise source are apparent. 
However, in the limit of Poissonian white noise (A' ^ A) such a distinction disappears; 
under such a condition equation fITOl) becomes 

CiT) = ^e-''\ (71) 
m 

Unless the relaxation of the coUisional effects is relevant, these effects and those caused 
by the surface should be indistinguishable and equation (I7T|) would describe accurately 
the loss of correlation, with f given in terms of rj rather than 7 and/or A separately. 
The expression for the intermediate scattering function resulting from equation (17T!) 

is 

J(Air, t) = exp [-x^ (e-^* + r^t - l)] , (72) 
where the so-called shape parameter x [301 EI] is defined as 

X' ^ {vDAK'/v'- (73) 



From this relation we can obtain both the mean free path I = f a/ (fg) and the self- 
diffusion coefficient D = fiv^) [Einstein's relation). When the coverage increases, the 
collisions among adsorbates also increase, and so A (see section 12.2.4^ and therefore 77. 
As can be easily shown, equation (1721) displays a Gaussian decay at short times that does 
not depend on the particular value of 17, while at longer times it decays exponentially 
with a rate given by ri~^. Thus, with rj, the decay of the intermediate scattering function 
becomes slower. The type of decay is important concerning the width and the shape of 
the dynamic structure factor [T2l 130]. 
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The above described effects can be quantitatively understood by means of the 
expression of the dynamic structure factor obtained analytically from equation fl72l) . 



Here, the F symbol in the first line denotes both the Gamma and incomplete Gamma 
functions (depending on the corresponding argument), respectively. As can be noted in 
the high friction limit, equation (17^ becomes a Lorentzian function, its full width at half 
maximum (FWHM) being F = 2?7x^, which approaches zero as t] increases (narrowing 
effect). This in sharp contrast to what one could expect — as the frequency between 
successive collisions increases one would expect that the line shape gets broader (effect 
of the pressure in the spectral lines of gases). The physical reason for this effect could be 
explained as follows. As rj increases the particle's mean free path decrease and therefore 
correlations are lost more slowly. In the limit case where friction is such that the particle 
remains in the same place, the van Hove function becomes a 5-function, the intermediate 
scattering function remains equal to 1 and the dynamic structure factor consists of a 
5-function at = 0. Conversely, in the low friction limit the line shape is given by a 
Gaussian function, whose width is F = 2^2 In 2 A/fesT/m AK, which does not depend 
on rj. This is the case for a two-dimensional free gas [SUl [32] • This gradual change 
of the line shapes as a function of the friction and/or the parallel momentum transfer 
leading to a change of the shape parameter x is known as the motional narrowing effect 
pT| [T2l [30] . Note that in our approach the friction is related to the coverage. Thus, at 
higher coverages a narrowing effect is predicted for a fiat surface [T5j . 

2.2.7. The harmonic oscillator. Bound adparticles. In contrast with the case of a 
dynamics where V does not play a relevant role, we can devise a particle fully trapped 
within a potential well. The harmonic oscillator is an appropriate working model to 
understand the physics associated with this problem. This model also allows us to 
understand the behaviour associated with the T mode, which comes precisely from 
the oscillating behaviour undergone by the particle when the diffusional motion is 
temporarily frustrated. 

For a harmonic oscillator, the behaviour of the adparticle becomes very apparent 
when looking at the corresponding velocity autocorrelation function, which reads [T2|[M] 
as 



S{AK, u) = — x"'^' Re jx-'^^/^lf (x' + lu h) - f(x' + lu/r,, x')]| 
_ e^' ^ (-1)"X^" (X^ + n)7] 




(74) 




(75) 



Here, 




(76) 
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uoq being the harmonic frequency. Alternatively, equation fl75p can also be recast as 

m u 

with 

5 = (tan)-i (^^^ . (78) 

Note that equation ( 17Ti) can be easily recovered after some algebra in the limit uq ^ 
from either equation (175!) or (1771) . In the case of anharmonic potentials, provided that 
we work within an approximate harmonic regime, ujq would represent the corresponding 
approximate harmonic frequency. 

The only information about the structure of the lattice is found in the shape 
parameter through AK [see equation fl751) ]. When large parallel momentum transfers 
are considered, both the periodicity and the structure of the surface have to be taken into 
account. Consequently, the shape parameter should be changed for different lattices. 
The simplest model including the periodicity of the surface is that developed by Chudley 
and Elliott |33j, who proposed a master equation for the pair-distribution function in 
space and time assuming instantaneous discrete jumps on a two-dimensional Bravais 
lattice. Very recently, a generalized shape parameter based on that model has been 
proposed to be [31] 



where, for our approach, it has been written t] instead of 7. Here, ri,(AK) represents 
the inverse of the correlation time and is expressed as 

r.(AK) = z/^Pj [1 - cos(j ■ AK)], (80) 
j 

u being the total jump rate out of an adsorption site and Pj the relative probability that 
a jump with a displacement vector j occurs. 

Introducing now equation fl77|) into ([9]) leads to the following expression for the 
intermediate scattering function 

I(AK, t) = exp <^ [cos 5 - e-''*/^ cosicot -5)]\. (81) 

The argument of this function displays an oscillatory behaviour around a certain value 
with the amplitude of the oscillations being exponentially damped. This translates 
into an also decreasing behaviour of the intermediate scattering function, which also 
displays oscillations around the asymptotic value. This means that after relaxation the 
intermediate scattering function has not fully decayed to zero unlike the free-potential 
case. Again, in the limit Uq ^ 0, equation (ISTI) approaches equation (1721) . 
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In order to obtain an analytical expression for the dynamic structure factor, it is 
convenient to express equation flHTj) as 

/ 1 \ m+n 



m! n\ 

m,n=0 



where 



with 



X g«('"-")'^g-(™+")'7*/2g«(ni-n)wi j^g2^ 



_ ^0 r/^{2(r//2)c^sin^ + - (r//2)^] cos^} 
' iU [(r//2)2 + c^2]2 ' l^^^ 

A3 = —T^ri^ — 85 

uj 2(?7/2 — tuj)^ 

- ^ 2(7^/2 + ^^^^ 

where the coefficients Ai has been put in terms of 77, a) and 6. From equation (!82l) . it 
is now straightforward to derive an expression for the dynamic scattering factor, which 
results 

71 ^-^ mini 

m,n=0 

X {m + n)r]/2 

[cu- (m-n)a;]2 + [(m + n)r//2]2' ^ ^ 



For a harmonic oscillator, there is no diffusion and, therefore, equation flHTI) is only valid 
when m ^ n. All the Lorentzian functions contributing to equation fl87|) are due to the 
creation and annihilation events of the T mode. These Lorentzians are characterized by 
a width given by F = {m+n)'r]/2, which increases with rj. This broadening (proportional 
to 1]) undergone by the dynamic structure factor is thus contrary to the narrowing effect 
observed in the case of a fiat surface. It can be assigned to the confined or bound 
motion displayed by the particle ensemble when trapped inside the potential wells. 
Hence, in order to detect broadening of the line shapes in surface diffusion experiments, 
adparticles must spend some time confined inside potential wells, since the broadening 
will be induced by the presence of temporary vibrational motions. 

2.2.8. General periodic surface potentials. Temporary trapped adparticles. As seen 
above, the broadening of the diffusion line shapes is provoked by the temporary 
trapping. To demonstrate this assertion here we are going to consider a general velocity 
autocorrelation function which keeps the functional form of equation (1771) . but whose 
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parameters do not hold the same relations as those characterizing a harmonic oscillator 
[12] . That is, 

C{t) = ^ e"'^" cos(cDt + 5), (88) 

where the values of the parameters fj, uj and 5 are obtained from a fitting to the numerical 
results — there is no any relation among them as in the case of the harmonic oscillator 
model. From equation (IHH!) one easily reaches the corresponding expression for the 
intermediate scattering function, 

J(Air, t) = e-^?^>'*) 



^-^ mini ' * 



m,n=0 



^ ^i{m-n)S ^-(m+n)f)t ^i(m-n)ut ^gg^ 

which is analogous to equation ( l82l) . In equation ( !89l) . 

/(cD, t) = ii + iat + ise'^e-^''-^^)* + i4e-^^'e-(^^+^^)*, (90) 

and 

~ f/^[2?7u;sin5 + (cj^ — r]^) cos5) , , 

^1 = ' ^ ^ 

~ _ ^^(??cos(^-tI;sin(5) 

~2 

is = 7^7^^^, (93) 

~2 

i4 = ^t-l-.r (94) 

Unlike the case of the harmonic oscillator, notice now that there is a linear dependence 
on time in / because of the independence of f], uj and 5. This leads to an exponentially 
decaying factor in equation (l89ll . which accounts for the diffusion and that makes 
the intermediate scattering function to vanish at asymptotic times. In this sense, 
the intermediate scattering function can be considered as containing both phenomena, 
diffusion and low vibrational motions. This effect is better appreciated in the dynamic 
structure factor, 

S(AK,uj) = - y L^^;^f(-+")i-i-e*(— 

vr m\n\ 

m,n=0 

^ XM2 + {m + n)f] ^g^^ 

[uj — {m — n)ujY + [x^^2 + ("^ + n)fiY 
This general expression clearly shows that both motions (diffusion and oscillatory) 
cannot be separated at all. The Q-peak is formed by contributions where m = n, 
for which each partial FWHM is given by 

TQ = XiA2 + 2mfi. (96) 
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Analogously, the T-peaks come from the sums with n ^ m and each partial FWHM is 
given by 



If the Gaussian approximation is good enough, the value of fj will not be too different 
from the nominal value of rj and, therefore, both peaks will display broadening as rj 
increases. This is a very remarkable result since a relatively simple model, as the one 
described here, can explain the corresponding experimental broadenings observed with 
coverage. Thus, broadening arises from the temporary confinement of the adparticles 
inside potential wells along their motion on the surface [15]. The problem of the 
experimental deconvolution has been already discussed elsewhere [31]. Here we would 
like only to mention that using this simple model, such deconvolutions would be more 
appropriate in order to extract useful information about diffusion constants and jump 
mechanisms. Finally, as mentioned before, the motional narrowing effect will govern the 
gradual change of the whole line shape as a function of the friction or, equivalently, the 
coverage, the parallel momentum transfer and the jump mechanism. 

3. Results 

For a better understanding of the concepts introduced in section [21 here we present 
results for two different types of two-dimensional surfaces: fiat and corrugated. 
Since much research has been developed for Na atoms adsorbed on a Cu(OOl) 
surface, we will carry out our numerical simulations taking into account the periodic, 
separable corrugated potential found in the literature for this system (see, for example, 
reference [7]). 

3.1. Computational details 

As in reference [7j, we have considered two values of the coverage in our calculations, 
6 = 0.028 and 0.18, where = 1 corresponds to one Na atom per Cu(OOl) surface atom 
— or, equivalently, a = 1.53 x 10^^ atom/cm^ [8]. All the friction coefficients are given 
in atomic units. Taking into account the values a = 2.557 A and p = 2 A for the unit 
cell length and the Na atomic radius, and using equation (139|) . the values of A associated 
with the coverages used here are A = 3.34 x 10~^ a.u. for 6 = 0.028 and A = 2.15 x 10~^ 
a.u. for = 0.18 at T = 200 K. For the surface friction we have considered the value 
also given in the literature [7|, 7 = O.Icjq = 2.2049 x 10^^ a.u., where uq is the harmonic 
frequency associated with the periodic surface potential (see section l373l) . With these 
two frictions (A and 7) the resulting total friction is rj ^ 2.53 x 10"^ a.u. for 6 = 0.028 
and 1] ^ 4.35 x 10~^ a.u. for 6 = 0.18. Regarding the collision relaxation rate, we have 
assumed A' = 10"^ a.u., which allows us to use the standard Langevin equation with a 
white shot noise. As far as we know, no information about that parameter is found in 
the literature and therefore the interacting single adsorbate approximation could also 




(97) 
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give us an estimation of the duration of the adsorbate-adsorbate colhsion. Finally, the 
Verlet algorithm [35] is used to solve the corresponding Langevin equation. 

3.2. Flat surface model. The low corrugation regime 

First we analyze the case of adsorbate diffusion on a flat surface, i.e., the case of a 
two-dimensional gas [301 132] • This example is representative of low corrugated surfaces, 
where the role of the activation barrier is negligible. 

As seen in section I2.2.4[ when = 0, the dynamical magnitude {x'^{t)) displays 
two different dynamical regimes: ballistic (t <^ l/r^) and diffusive {t ^ ^/v)- This can 
be seen in figure lU^a), where {x^{t)) is effectively proportional to t in the long time 
regime and depends on t^ at short times (of the order of 1/r/). In the linear regime, 
the slope of increases as rj decreases, i.e., in accordance with Einstein's law (see 

equation (|68|) ). the diffusion decreases with the friction. By fitting the linear part of 
the graphs (dotted lines) to equation ( l68i) . we get D = 6.045 x 10~^ a.u. for 9 = 0.028 
and D = 3.510 x 10~^ a.u. for 6 = 0.18. Note that these diffusion coefficients are not 
only related to the friction due to the surface (i.e., the adsorbate-substrate coupling), 
but also to the collisions among adsorbates. Indeed, we can see that for a given surface 
friction, the diffusion is inhibited by such collisions — it decreases with the coverage. 
This is a very remarkable result since in real experiments the surface friction is fixed 
and diffusion can then be studied taking only into account the coverage of the surface. 
On the other hand, from these values for D we obtain, respectively, rj = 2.50 x 10^^ 
a.u. and rj = 4.35 x 10~^ a.u., which are in a good agreement with those used in our 
calculations (see section [3TT1) . This agreement between the simulation and the analytical 
model is particularly important because it validates our standard Langevin description. 

The effect of the collisions also produces an important effect on the velocity 
autocorrelation function, C{t), plotted in Fig. |l](b). As seen, as 9 increases this function 
decays faster. In particular, from the fitting of these results to equation flTTl) (normalized 
to unity) we have obtained rj ^ 2.52 x 10~^ a.u. for 6 = 0.028 and rj ^ 4.32 x 10~^ 
a.u. for 6 = 0.18, which again show a good agreement with the actual values employed 
in our simulations. This fast decay also indicates that higher order correlations (e.g., 
correlations at three or four times) will decay much faster. Their effect will be then 
negligible on the intermediate scattering function, thus validating the use of the Gaussian 
approximation, since such higher order correlations will not be relevant when passing 
from equation ([8]) to ([9]). The validity of the Gaussian approximation can be seen 
more explicitly when comparing the intermediate scattering functions obtained from the 
calculations with that fitted using equation fl72|) . These results are plotted in figure IH^c) 
for AK = 1.23 A^^, where the fitting has been carried out with the values of rj used 
in the simulation and considering x as the fitting parameter. We observe that not only 
the correspondence between the simulations and the analytical formulas are excellent 
from the figure, but also from the fitted values of x- Xfu = 3.16 vs Xsim = 3.16 for 
6 = 0.028, and Xfu = 1-81 vs Xsim = 1-83 for 9 = 0.18. Notice that in agreement with 
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Figure 4. Dynamical magnitudes for two different values of the coverage, 9 = 0.028 
(blue) and 9 = 0.18 (red): (a) (^^(i)), (b) C{t), (c) /(AK,i) and (d) S'(w), with 
Ai^ = 1.23 A-i. Dotted lines are the numerical fittings to the corresponding analytical 
formulas given in section [2l 



equation (!72l) . I{t) displays an initial Gaussian falloff at short times, while for longer 
times its decay is exponential. Moreover, as also expected from equation fl72l) . and 
equation (1731) . as the coverage increases a slower decay is observed. 

Finally, in figure 111(d) we have plotted the dynamic structure factor after time 
Fourier transforming the intermediate scattering function obtained from both the 
numerical calculations and their corresponding analytical fittings. Note that, a 
narrowing of the Q peak is predicted with the coverage [15]. Regarding the shape 
of this peak, it can be shown that it is a mixture of Gaussian and Lorentzian functions. 
The Gaussian behaviour is ruled by the short time limit of the intermediate scattering 
function, while the Lorentzian behaviour arises from the long time exponential decay. 
The shape of the Q peak will depend on which of the two extreme regimes is dominant 
(motional narrowing effect) This analysis should be carried out when experimental 
results are deconvoluted [30] since it is very common to see fittings of the Q peak to a 
pure Lorentzian function. 

3.3. Periodic, corrugated surface model 

Now we are going to study the case of a periodic surface model described by the potential 

V{x, y) = Vq[2- cos(27r2;/a) - cos(27ry/a)] , (98) 

with a = 2.557 A and Vq = 37.32 meV. In this case, the surface corrugation is relatively 
strong and cannot be neglected, thus inducing important effects in the adsorbate 
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Figure 5. Velocity autocorrelation function, C{t), for two different values of the 
coverage: (a) = 0.028 and (b) = 0.18. Solid lines denote the numerical results 
obtained from the simulations and dotted lines the fittings to equation ([55]) . 



dynamics. These effects can be well understood by using the third model described 
in section I2.2.8[ 

After figures [5l the diffusive and trapped regimes that appear under the presence 
of the surface corrugation give rise to the oscillating behaviour of the velocity 
autocorrelation function. The adparticles trapped inside the potential wells give rise 
to the oscillations, while the exponential damping arises as a consequence of the loss of 
correlation with time of both the particles moving outside the wells and those trapped 
inside. This behaviour thus approaches that described by equation (188|) . as seen from 
the fittings (dotted line) with parameters fj = 5.13 x 10~^ a.u., u = 2.14 x 10~^ 
a.u. and 6 = -0.874 for 6 = 0.028 (see figure E^a)), and fj = 5.96 x 10"^ a.u., 
Cu = 2.05 X 10-^ a.u. and S = -0.776 for 9 = 0.18 (see figure El^b)). Note that, 
since we have contributions from running and trapped or bound trajectories, here the 
fitting is not so good as that seen in section [3^ for V = 0. Recently, it has been shown 
[15] that if the velocity autocorrelation function is expressed as a linear combination 
of the corresponding functions for free particle and anharmonic oscillator, the fitting 
is highly improved. Obviously, the Gaussian approximation is no longer valid due to 
the presence of the corrugated potential. Nonetheless, it is important to stress that the 
model given by equation (|88ll is still valid to interpret the results obtained; not only the 
fitted parameters are of the order of those employed in the calculations, but they also 
follow the trend that one might expect from the exact ones. That is, as fj increases the 
exponential damping in C{t) is stronger and a slight shift of the oscillations is observed. 

As seen in figure |6l the effects of the trapping or bound motion are also apparent 
in the intermediate scattering function, which is plotted for AK = 1.23 A~^, along two 
different directions and two coverages: panels (a) and (a') for 6 = 0.028 along the [110] 
and [100] directions, respectively; and panels (b) and (b') for 6 = 0.18 and along the 
same two directions. Note first that again here we can distinguish two decay regimes 
(which are present along both directions). The long time regime in the upper panels 
is characterized by an exponential damping, as in the case for V = 0, which is typical 
of diffusion and will give rise to a predominance of the Lorentzian-like behaviour in 
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Figure 6. Intermediate scattering function for AX = 1.23 and two different 
values of the coverage: 9 = 0.028 in panels (a) and (a') along two directions, [110] and 
[100] and = 0.18 in panels (b) and (b') along the same two directions. In the bottom 
panels observe details of the short time regime. Solid lines denote the numerical results 
obtained from the simulations and dotted lines the fittings to equation ([55]) . 

the dynamic structure factor (see below). On the other hand, the short time regime 
in the bottom panels clearly shows the damped oscillating behaviour of C{t). These 
oscillations at short times will give rise to the T-mode peaks associated with the low 
frequency vibrational motions induced by the potential wells in the trapped trajectories. 
Moreover, the oscillating behaviour is more patent with the coverage giving rise to an 
inversion in the behaviour of the long time tail of I{t) . Unlike the intermediate scattering 
function associated with a flat surface dynamics, this faster decay with coverage will 
lead to a broadening of the Q peak in the dynamic structure factor. Moreover, also 
notice that the timescales ruling the decay are about two orders of magnitude larger 
than in the case with V = This is due to the presence of the potential, which allows 
to keep the phase correlation in equation (IHl) for longer times due to the trapping inside 
the potential wells [T5] . 

Regarding the dynamic structure factor, displayed in figure [7] (as before, the four 
panels correspond to the same coverages and directions as in figure [6]), the presence of 
the potential has three important differences with respect to the fiat surface case. First, 
the trapping trajectories due to the surface corrugation give rise to the appearance of two 
symmetric peaks around the Q peak, the T-mode peaks which are located at —5 meV 
and 5 meV, approximately (in the lower panels of figure [7] only the right peak is showed); 
these peaks broadens with coverage. Second, the faster decay of I{t) as the coverage 
increases translates into a broadening of the Q peak, which is in accordance to what one 
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Figure 7. Dynamic structure factor for Aif = 1.23 and two different values of 
the coverage: 9 = 0.028 in panels (a) and (a') along two directions, [110] and [100] 
and 9 = 0.18 in panels (b) and (b') along the same two directions. In the bottom 
panels observe details of the peaks corresponding to the T mode. Solid lines denote 
the numerical results obtained from the simulations and dotted lines the time Fourier 
transform of the fittings to equation ([89]) . 

can observe experimentally [?]. Third, the Q peak is narrower than the corresponding 
in the fiat surface case, around two orders of magnitude. This narrowing is due to the 
inhibition of the diffusive motion induced by the corrugation of the potential (a large 
number of particles get trapped inside the wells). In any case, the diffusion coefficient 
is governed by the Einstein relation and therefore decreases with the coverage. 

Finally, it is worth stressing that there is a strong correlation between diffusion 
and low frequency vibrational motions [I2]. This fact is easily appreciated in the two 
bottom panels of figure [71 where a strong overlapping of the two peaks is noticeable. In 
the Gaussian approximation, this overlapping can be understood from equation (|95l) . 
In general, in order to extract the width of the Q and T peaks, experimentalists usually 
follows a deconvolution procedure where an effective Lorentzian function is assumed for 
each peak independently. After the study presented here (as well as in references [121 El] 
within the single adsorbate approximation), this deconvolution procedure could lead to 
inconsistencies when extracting information because both peaks are strongly overlapped. 
Indeed, it can be shown [36j that depending on the parallel momentum transfer AK 
considered, the overlapping of both peaks may become so strong that the Lorentzian 
fitting procedure would not be valid anymore. In the light of this study we propose a 
more general fitting procedure based on equation (!95l) . 
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4. Conclusions 

In this work, we liave presented a full stochastic description of surface diffusion and low 
frequency vibrational motion for adsorbate/substrate systems with increasing coverages. 
A Langevin equation is used to describe the adparticle dynamics, where a Gaussian 
white noise simulates the effects of adsorbate-substrate friction and a white shot noise 
is used for the adsorbate-adsorbate interactions within what we call the interacting 
single adsorbate approximation. As shown, this theoretical stochastic framework not 
only provides a simple and complementary view of the processes taking place on a 
surface, but also explains the experimental observations, i.e., the line shape broadening 
of the Q peak ruling the diffusion process as a function of the coverage. The idea of 
replacing the dipole-dipole interaction by a white shot noise turns out to be crucial 
because for long time processes, as the diffusion one, a high number of collisions occurs 
and the statistical limit seems to wipe out any trace of the true interaction potential. 

The main points of this work can be emphasized as follows. First, the simplicity of 
the stochastic approach developed here is in sharp contrast to the Langevin molecular 
dynamics simulations commonly used to study adsorbate diffusion. This treatment 
allows to derive analytical formulas that render some light into the physics involved in 
surface diffusion, something that is not possible by means of massive molecular dynamics 
calculations which provide numbers that cannot be easily handled by using an analytical 
data processing methodology. Second, by using this numerical and analytical treatment, 
we have been able to show that the main reason for the broadening of the Q peak is 
the trapping induced by the surface corrugation. Third, this analysis constitutes a 
robust and simple methodology that can be of great help for the experimentalists in the 
interpretation of their results. And fourth, the results obtained here can be of relevance 
in problems where friction has to be varied gradually. Note that the surface friction 
is a fixed parameter, since it depends only on the type of surface. However, relating 
the coverage with a friction (A) seems to be an appropriate manner of studying such 
variations. 

Finally, this class of studies can obviously not replace further investigation at 
microscopic level and calculations from first principles of adsorbate dynamics. For 
example, the Lau-Kohn long range interaction via intrinsic surface states [37] has 
been proposed to use a coverage dependent surface electronic density modifying the 
corrugated potential. It has also been used to explain the remarkable increasing of the 
T-mode frequency with coverage and where the dipole-dipole interaction is not able to 
reproduce this behaviour. Our simple stochastic model provides a complementary view 
of the diffusion and low-frequency vibrational motion features observed as peaks around 
or near zero energy transfers (or long time regime) in the scattering law. 
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